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In a recent Letter entitled "A new numerical method for obtaining gluon distribution functions 
G{x,Q^) = xg{x,Q^), from the proton structure function F2^{x,Q'^)" we derived an accurate 
and fast algorithm for numerically inverting Laplace transforms, which we used in obtaining gluon 
distributions from the proton structure function F2^{x, Q^). We inverted the function g{s), where s 
is the variable in Laplace space, to G{v), where v is the variable in ordinary space. Since publication, 
we have discovered that the algorithm does not work if g{s) — s- less rapidly than 1/s, as s — s- oo. 
Although we require that g{s) — >■ as s -t- oo, it can approach as with < /3 < 1, and still be 
a proper Laplace transform. In this note, we derive a new numerical algorithm for just such cases, 
and test it for g{s) — the Laplace transform of -j^. 



I. INTRODUCTION 

In an earlier note [l[ we developed an algorithm to numerically invert Laplace transforms, in order to find an analytic 
solution for gluon distributions. We used a global parameterization of the proton structure function, F2^{x, Q^) and 
a LO (leading-order) evolution equation for . Subsequently, going to NLO (ncxt-to- leading order) in the strong 
coupling constant as{Q^), we have discovered the algorithm failed badly. Detailed investigation showed that the cause 
of the problem was that for this g{s)-whose Laplace transform was our desired NLO gluon distribution G{v) — it went 
to slower than 1/s as s — oo, where s is the Laplace space variable. The purpose of this note is to derive a new 
algorithm for such cases, not covered in [l[, which can be modeled by Laplace transforms of the type 

g{s) ^, < /? < 1. (1) 

To illustrate the new method, we will numerically invert g{s) = the Laplace transform of and test its accuracy, 
as well as show the inadequacy of the original algorithm for such a case. 

II. NUMERICAL INVERSION OF LAPLACE TRANSFORMS 

Let 17(5) be the Laplace transform of G{v). The Bromwich inversion formula for G{v) is given by 

-1 f c-\-i 00 

G{v)^C-'[g{s)-v] = — g{s)e^'ds, (2) 

where the real constant c is to the right of all singularities of g{s). Making an appropriate coordinate translation in 
s so that c = 0, we can write 

1 /• +ioo 

G{v)^£-'[g{s);v] = — g{s)e"Us. (3) 

Our goal is to numerically solve Eq. ([3]). The inverse Laplace transform is essentially determined by the behavior of 
17(5) near its singularities, and thus is an ill-conditioned or ill-posed numerical problem. We suggest in this note a 
new algorithm that takes advantage of very fast, arbitrarily high precision complex number arithmetic that is possible 
today in programs like Mathematica , making the inversion problem numerically tractable. 
First, we introduce a new complex variable z = vs and rewrite Eq. ^ as 

G{v) = —- / g[-)e^dz. (4) 



In our earlier paper we proceeded to make a rational approximation to the exponential e^, under the tacit 
assumption that g{s) went to sufficiently rapidly as s went to 00. This is not a valid assumption if g[s) goes to as 
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l/x'^, < /? < 1. In this case, we will rewrite Eq. ^ as 
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oo 




|^).Vd. (5) 



and make a rational approximation to z e^, using the Pade approximant 



Pade[ze^O, (2Ar-3,2iV)], (6) 

where the numerator of the function Pade in Eq. © is a polynomial in z of order 2A'' — 3 and whose denominator is 
a polynomial in z of order 2iV, and the expansion is around z = 0. Let ai, i = 1, 2, . . . , 2N be the 2N complex roots 
of the denominator, i.e., its 2N complex poles, and let uji be the complex residues of the Pade. It can be shown that 
they have the following properties: 

1. Re ai > 0, so that its poles are all in the right-hand half of the complex plane. 

2. A'^ distinct complex conjugate pairs of the complex numbers {uti, ai), such that the sum of the fc"^ pair, 



z- ak z ~ ak 

is real for all real z. 



(7) 



3. The integrand vanishes faster than 1/ R as i? — > oo on the semi-circle of radius R that encloses the right hand 
half of the complex plane, since the approximation vanishes as and g{s) that corresponds to a non-singular 
G{v) also vanishes for R —J- oo. 

Since g(z/v) must vanish for z — > oo and our approximation for z^e^ in Eq. (|6]) vanishes for z — > oo, we can form a 
closed contour C by completing our integration path of the modified Bromwich integral in Eq. ^ with an infinite 
half circle to the right half of the complex plane. As mentioned earlier, g(z/v) has no singularities in this half of the 
complex-Z plane. It is important to note that this contour is a clockwise path around the poles of Eq. ([6]), which come 
from our approximation to ze^. What we need is the negative of it, i.e., the contour —C which is counterclockwise, 
so that the poles are to our left as we traverse the contour —C. Therefore, we rewrite Eq. (g]) as 

G{v) « ^— / ^-^Pade[zV, 0, (2M - 3, 2M)] 
2Triv J(j z^ 



/ ^^Pade[z2e",0, (2M-3,2M)] 

(8) 
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To obtain Eq. ([5]), the final approximation to G{v), we used Cauchy's theorem to equate the closed contour integral 
around the path —C to 2-Ki times the sum of the (complex) residues of the poles. Since the contour —C restricts us 
to the right-hand half of the complex plane, no poles of g(z/v) were enclosed, but only the 2N poles ai of the Pade 
approximant of z^e^. Using the properties cited above of the complex conjugate pairs — {uii,ai) and (oji,ai) — after 
taking only their real part and multiplying by 2-we have simultaneously insured that G{v) is real , yet only have had 
to sum over half of the residues. The above scheme now works for the slowly converging case of g{s) = l/s*^, where 
< /3 < 1, because we are dividing it by z^ to make the quantity that we are inverting , i.e., g{z/v)l z, go to more 
rapidly than 1/s as s — )• 0. 

Two properties of Eq. ([8]) arc worth emphasizing here: 

1. The AN coefficients [ai^uJi) are complex constants that are independent of v, only depending on the value of 
2N used for the approximation, so that for a given 2N , they need to be evaluated only once — in essence, they 
can be tabulated and stored for later use. 

2. The real parts of the residues uJi alternate in sign and are exceedingly large — even for relatively modest 2N , 
making round-off a potentially serious problem. Thus, exceedingly high precision complex arithmetic is called 
for, often requiring 60 or more digits. However, this is not a serious problem — either in speed or complexity of 
execution — for an algorithm written in a program such as Mathematica • 

A concise inversion algorithm in Mathematica that implements Eq. (jS]) is given in Appendix [XI 
We will now test the accuracy of our numerical Laplace inversion algorithm by using g(s) = the Laplace 
transform of G{v) = ^/v. 
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III. COMPARISON OF EXACT SOLUTION AND NUMERICAL LAPLACE INVERSION RESULTS 

To illustrate our inversion routine for slowly converging Laplace transforms, we plot in Fig. [l]the numerical inversion 
of = the Laplace transform of G{v) = The solid curve is the function G{v) = and the (blue) dots 
are a result of using the algorithm in Appendix \^ using Mathematica 0, with 2N 20. As seen, the agreement 
is excellent. Also shown in Fig. [1] are the (red) squares, which arc the result of our using our earlier algorithm ll| 
for the calculation of the inverse Laplace transform. Clearly, this is a disaster. As detailed in our earlier work 
the original algorithm was exact if the function in v space was a polynomial of degree < AN — 1, so that it is a very 
powerful tool when used properly, i.e., with a g{s) that goes to as fast or faster than 1/s as s — > oo. However, for 
5'(s)'s that don't vanish this rapidly, that algorithm is deficient and must be replaced by the new algorithm given in 
this note. 

As a practical matter, if we have a g{s) to invert that is purely numerical and whose behavior at oo is completely 
unknown (except for the requirement that it mitst vanish at oo), we must test its convergence properties by trying both 
algorithms. If they agree reasonably numerically, we should continue using the more powerful algorithm of Ref. 
It they disagree completely, as is the case for the squares in Fig. [U then one should use the new algorithm developed 
in this note. 

We will discuss in Appendix |X] the inherent accuracy of the algorithm used to get the dots in Fig. [T] As noted in 




FIG. 1: A plot of G{v) = vs. V. The solid curve is G{v). The (blue) dots are the approximate values calculated from 

Eq. ([HI , using 2N =20 and g(s) = the (red) squares are the result of using 2A'^ = 20 in the algorithm described in Ref. 
which should only be used for Laplace transforms that go to at oo as fast, or faster than 1/s. 

[l[ , our inversion routine clearly has a wide variety of additional applications in solving both integral and differential 
equations. 



IV. CONCLUSIONS 



We have achieved arbitrarily high accuracy in numerically inverting Laplace transforms of the type g{s) = < 
/3 < 1. As an example, we show in Fig. [T]a comparison of the solution G{v) — i.e., the inverse Laplace transform 

of g{s) = with the exact answer, and demonstrate the algorithm's inherent accuracy in Fig. [21 showing that it is a 

power law in the reciprocal of the number of terms used in the expansion, and thus, of arbitrary accuracy. Because 
the assumed g{s) go to very slowly as s 00, the algorithm that we previously developed in Ref. [l| does not work 
here, and we discuss when and when not to use it. We are currently using our new algorithm (this communication) 
to analytic decoupled solutions for LO singlet structure functions and gluons Q , as well as for obtaining NLO gluon 
distributions. 
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Appendix A: Mathematica Laplace Inversion Algorithm 

NlnverseLaplaceTransf ormBlock2 [g_, s_, v_, twoN_,prec J : =Module [ 
{ Omega, Alpha, M,p, den, r,num} , 

(M=2*Ceiling[twoN/2] ;If [M<6,M=6] ;p=PadeApproximant [ z"2*Exp[z] ,{z,0,{M-3 ,M}}] ; 
den=Denominator [p] ;r=Roots [den==0,z] ; Alpha=Table [r [ [i ,2] ] ,{i,l,M}] ; 
num=Numerator [p] ;hospital=nmn/D[den,z] ; 

Oinega=SetPrecision [Table [hospital/. z->Alpha[[i]] ,{i,l,M,2}] ,prec+50] ; 
Alpha=SetPrecisioii [Table [Alpha[[i]] ,{i,l,M,2}] ,prec]) ; 

SetPrecision[-(2/v) Sum [Re [Omega [[i]] (g/ (v*s) ~2) / . s->Alpha [ [i] ] /v] ,{i,l,M/2}] ,prec] 

] 

In the above algorithm, g = .9(s), s ~ s, v ~ v, t'wo'N=2N in Eq. ([S]), and prec — desired precision of calculation. 
Typical values are twoN = 10 and prec = 70. The algorithm, which is quite fast, returns the numerical value of G{v). 

The algorithm first insures that M=twoN is even and > 6. It then constructs p, the Pade approximant whose 
numerator is a polynomial in z of order M — 3 and denominator a polynomial in z of order M, whose expansion 
is around z ~ 0. It finds r, the complex roots of the denominator, which are a^, the poles of Eq. Using 
L'Hospital's rule, it finds the residue uji corresponding to the pole a^. At this point, all of the mathematics is 
symbolic. It next finds every other pair of (0^,0;^) to the desired numerical accuracy; they come consecutively, i.e., 
Qfi = (52, ^i^i = <^2, cea ~ oii, "^3 = <^4, etc. Finally, it takes the necessary sums, again to the desired numerical 
accuracy, but only over half of the interval i = 1, 3, . . . , A^, by taking only the real part and multiplying by 2. 
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FIG. 2: A log-log plot of the fractional accuracy of the numerical Laplace inversion of the original function G{v) = vs. 
twoN, with twoN=2A'^, the order of the rational approximation in Eq. ((HJ. 

If g{s), the input to the algorithm, is an analytic relation and w is a pure number (from the point of view of 
Mathematica, 31/10 is a pure number, but 3.1 is not), then, for sufficiently high values of prec, you can achieve 
arbitrarily high accuracy. If we define the fractional accuracy as 1 — G(f )numGricai/G'(w)truoj numerical tests on many 
different functions shows that it goes to for large 2N as a power law in 1/2N. We illustrate this in Fig. [51 where 
we show a log-log plot of the fractional accuracy vs. 2N. The straight line shows that the fractional accuracy is a 
power law in 2N, allowing one to achieve arbitrary accuracy in the inversion of g{s) = Further discussion of 
applications and numerical properties can be found in Ref. [l|. 
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